%% instructions: 

% to replicate baseline results: just run the file

% Notice: in absence of the MatLab parallel toolbox, please replace the
% parfor loop in VARirtrue_sim1shock_histories_futureN.m with a for loop.

%%
clear;  clc; close all;
rand( 'seed',1234)
randn('seed',1234)

addpath('../LibForSerial/LibCaggianoEtAl/','../LibForSerial/LibCaggianoEtAl/I-VAR tbx_/','../LibForSerial/LibCaggianoEtAl/I-VAR tbx_/utilsCesaBianchi/','../LibForSerial/LibCaggianoEtAl/I-VAR tbx_/utils4Figures/','../LibForSerial/LibraryLeSage/')

fast=0; %1 to be fast (avoid extraction future shocks in obtaining the GIRF - is a fairly good approximation/ 0 to consider them
% Notice that the fast option has to do only with confidence intervals (point estimates anyway are the proper ones)
% Notice: results in the published version have all be obtained with the fast=0 option
Nconsec = 2; % # of consecutive shocks, run it for 2,3,4 
SIMPLEversion = 0;      

%% COMPATIBILITY WITH Matlab versions
comp=1;  %put 0 if Matlab R2012b or R2014b. For alternative versions just try with 0 or 1, even though it is not said to work.
         %put 1 for MAC
         
%% Baseline data. Sample 1985M1 - 2019M12
  frequency = 12; % Monthly
load BBDdata

clearvars -except logEPU ...
                  fast comp projectdirDB projectdir Consec SIMPLEversion UncOrderedFirst Nconsec 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sample=[ 1985 1; 2019  12]; 
 
data    = xlsread('../Dataset/VARdataBBupdatedMonthly.xlsx','MainData');
dataInd = xlsread('../Dataset/VARdataBBupdatedMonthly.xlsx','IndicatorsContinuous');
startsample=find(data(:,1)==sample(1,1)*100+sample(1,2));
endsample  =find(data(:,1)==sample(2,1)*100+sample(2,2));

data = [data [NaN(12,1); data(13:end,4)./data(1:end-12,4)]];

data    = data(   startsample:endsample,2:end);
dataInd = dataInd(startsample:endsample,2:end);
clear startsample endsample 

 ADS          = dataInd(:,1);
 CFNAI        = dataInd(:,2); % short history, starting only in 1967
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA TRANSFORMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data(:,[3:11 12]) = log(data(:,[3:11 12]))*100;
%             LevelVsFirstDiff='IPGLevel' ;
              LevelVsFirstDiff='IPGGrowth';
if     strcmp(LevelVsFirstDiff,'IPGLevel')
    posIP = 3;
elseif strcmp(LevelVsFirstDiff,'IPGGrowth')
    posIP = 12;
end
%--------------------------------------%
dataMatrix=     [logEPU, data(:,1)          data(:,posIP)             data(:,10)]; VARspec = 2;

  data=dataMatrix;    
 clear dataMatrix
%-------------------------------------------------------------------------%
maxlag=6;
 %ConditioningVar = ADS;           stateTYPE = 1; THR = 0.0; %THR = max(ADS)  ;
  ConditioningVar = CFNAI;         stateTYPE = 1; THR = 0.0; %THR = max(CFNAI); 

%% specification
vardata = [data ConditioningVar lagmatrix(ConditioningVar.*data(:,1),[1:maxlag])];
labels={'UNC','Nominal Short-rate','Industrial production','Employment','State'};

 nlags=2; % Bootstrap does not work with nlags > 6 ... 
%THIS MUST MATCH nlags!!!
if     nlags==1
    labels_ex = {'Unc(-1)*StDe(-1)'}; 
elseif nlags==2
    labels_ex = {'Unc(-1)*StDe(-1)','Unc(-2)*StDe(-2)'};
elseif nlags==3
    labels_ex = {'Unc(-1)*StDe(-1)','Unc(-2)*StDe(-2)','Unc(-3)*StDe(-3)'};
elseif nlags==4
    labels_ex = {'Unc(-1)*StDe(-1)','Unc(-2)*StDe(-2)','Unc(-3)*StDe(-3)','Unc(-4)*StDe(-4)'}; 
elseif nlags==5
    labels_ex = {'Unc(-1)*StDe(-1)','Unc(-2)*StDe(-2)','Unc(-3)*StDe(-3)','Unc(-4)*StDe(-4)','Unc(-5)*StDe(-5)'}; 
elseif nlags==6    
    labels_ex = {'Unc(-1)*StDe(-1)','Unc(-2)*StDe(-2)','Unc(-3)*StDe(-3)','Unc(-4)*StDe(-4)','Unc(-5)*StDe(-5)','Unc(-6)*StDe(-6)'}; 
end
 
 
numend=size(labels,2); %number of the endogenous time series to take from the .xls file: remember to number them starting from the left

startdata=1; %1 to start from 1985:M1 

data_end   =vardata(startdata:end,1:numend);                % endogenous data
data_ex    =vardata(startdata:end,  numend+1:numend+nlags); % exogenous data + interaction terms (remember to put interaction variables last ordered in data_ex)
data_exFull=vardata(startdata:end,  numend+1:numend+maxlag);


[nobs] = size(data_end,1);

c_case=1;% 1 for constant, 2 for constant and trend in our VAR 

%/////////////////////////////////////////////////////////////////////////%
%to compute IRFs:
nsteps=36; 
    % NOTICE: in iterating the nonlinear VAR then we do that for "nsteps+5" times, before to save the ones needed. 
    % This can be useful to discard histories or sequences of future shocks that have started to explode (hence influencing the state-dep. girfs)
    % but that still cannot be recognized as such. However, in our case it is not important. 


pick=find(strcmp(labels,'UNC'));  % variable to be shocked (order as in the Cholesky decomposition)

ordend =pick;   % ordering (as in the Chol-decomp) of the I  endogenous variable with which to do the interactions (the shocked variable)
%#CHECK     -->has to be equal to pick
ordend2=numend; % ordering (as in the Chol-decomp) of the II endogenous variable with which to do the interactions (usually this is the same variable through which states are then defined)
%/////////////////////////////////////////////////////////////////////////%


% selection lags:
[laglength, AIC, logL] = InteractedVARlag(data_end,maxlag,c_case,data_exFull,0) % the same number of obs is maintained, the interactions change with the lag used

%#CHECK --> data_exFull = [];  % FOR LINEAR VAR
% OLS estimation Interacted-VAR
VAR   = VARmodel(data_end,nlags,c_case,data_ex);
VARprint(VAR,labels,labels_ex);

%-------------------------------------------------------------------------%
Sigma    = VAR.sigma;
[Cho, chol_flag] = chol(Sigma);
    S = Cho';
    if chol_flag~=0
        disp('-------------------------------------');
        disp('Error: VCV is not positive definite');
    end

residStruct=(S\VAR.residuals')';
%#CHECK
% figure(111)
% autocorr(residStruct(:,pick))
clear Cho chol_flag S Sigma

%% estimation nested linear VAR and nonlinear test
VARlin= VARmodel(data_end,nlags,c_case);

 LMstat=(nobs-nlags)*(log(det(VARlin.sigma))-log(det(VAR.sigma))) % see Canova ch. 4 
 dof=nlags*VAR.neqs ; % number of restrictions (additional interacted terms)
 pvalue=1-chi2cdf(LMstat,dof)

%%  algoritm to select the initial histories defining the two regimes
percL=10;
percH=90;
toler=5;  %tolerance threshold
    jj=0; %starting value for jj (Low  histories)
    zz=0; %starting value for jj (High histories)
    
%-------------------------------------------------------------------------%
if stateTYPE == 1
% 
%algoritm to select all histories for low and high uncertainty
%using tolerance regions, find threshold
     thrL1=prctile(data_end(:,ordend2),percL-toler);
     thrL2=prctile(data_end(:,ordend2),percL+toler);
     thrH1=prctile(data_end(:,ordend2),percH-toler);
     thrH2=prctile(data_end(:,ordend2),percH+toler);

%selecting histories:
LevAndChg = 1; % 1: Change < 0, Level < THR
               % 0:             Level < THR
observd=[0; diff(VAR.X(:,c_case+ordend2))];
for ii=1:VAR.nobs
    observ =     VAR.X(ii,c_case+ordend2);
if ~LevAndChg
% /////////////////             Level < 0 = Contraction ///////////////// %                     
    if observ<=THR  %
        jj=jj+1;
        historiesL(jj,:)=VAR.X(ii,:);
    else
        zz=zz+1;
        historiesH(zz,:)=VAR.X(ii,:);
    end
else
% ///////////////// Change < 0, Level < 0 = Contraction ///////////////// %   
    if observ<=THR  %
        if observd(ii) < 0
        jj=jj+1;
        historiesL(jj,:)=VAR.X(ii,:);
        end
    else
        zz=zz+1;
        historiesH(zz,:)=VAR.X(ii,:);
    end
end
end
%-------------------------------------------------------------------------%
elseif stateTYPE == 2
%-------------------------------------------------------------------------%
% 
%algoritm to select all histories for low and high uncertainty
%using tolerance regions, find threshold
     thrL2=prctile(data_end(:,ordend2),THR);
     
%selecting histories:
for ii=1:VAR.nobs 
    observ=VAR.X(ii,c_case+ordend2);
        
    if observ<=thrL2 %  <15th perc
        jj=jj+1;
        historiesL(   jj,:)=VAR.X(ii  ,:);
    else
        zz=zz+1;
        historiesH(   zz,:)=VAR.X(ii  ,:);
    end
end   
end
%#CHECK
% plot([historiesL(:,VAR.neqs+VAR.c_case)])
% plot([historiesH(:,VAR.neqs+VAR.c_case)])
if exist('historiesH','var') & exist('historiesL','var')
ToCheck = [quantile([historiesL(:,VAR.neqs+VAR.c_case)],[25 50 75]/100)
           quantile([historiesH(:,VAR.neqs+VAR.c_case)],[25 50 75]/100)];

       
    nhistL=size(historiesL,1);
    nhistH=size(historiesH,1); 
end
clear ii jj zz

%% obtain point estimates for state-dependent GIRFs 
%#CHECK --> ordend2 = 0; % FOR LINEAR VAR 
draws=500;
 typeImpulse=0; % 0 if 1st.dev shock, 1 if 1-unit increase shock. Whatever else gives exactly the size of the shock in st.dev. terms
 mode=0;        % 0 to draw from the empirical distribution of residuals for future shocks/ 1 from a Gaussian distribution respecting the VCV matrix for the future shocks

[OIRFavg1B,OIRF1B_optavg]=VARirtrue_sim1shock_histories_futureN(     VAR,nsteps+5,pick,historiesL,ordend,ordend2,draws,mode,typeImpulse,fast);
[OIRFavg2B,OIRF2B_optavg]=VARirtrue_sim1shock_histories_futureNcons( VAR,nsteps+5,pick,historiesL,ordend,ordend2,draws,mode,typeImpulse,fast,Nconsec);

% ---------- keep just the nsteps ahead needed
OIRFavg1B=OIRFavg1B(1:nsteps,:,:);
OIRFavg2B=OIRFavg2B(1:nsteps,:,:);
% ----------

% diagnostic explosiveness/instability:
disp('DIAGNOSTIC FOR GIRFS POINT ESTIMATES:')
disp(['- The number of explosive initial histories discarded is: ' num2str(OIRF1B_optavg.nexplos) ' for the Bad state w/ 1 consec shock '])
disp(['- The number of explosive initial histories discarded is: ' num2str(OIRF2B_optavg.nexplos) ' for the Bad state w/ 2 consec shocks'])
if max(OIRF1B_optavg.countRep)>0
disp('- It was needed to repeat some particularly extreme residual extractions for at least one initial history. For details see OIRF1_optavg.countRep')
end
if max(OIRF2B_optavg.countRep)>0
disp('- It was needed to repeat some particularly extreme residual extractions for at least one initial history. For details see OIRF2_optavg.countRep')
end
disp('------------------------------------------------------')

%% Obtain bootstrapped confidence intervals
perc=50; % percentile used to separe the two states inside the bootstrap procedure (all the FFR observations below it will be assigned to the low FFR (i.e. ZLB) state)
ndraws = 1000;
   %[INFavgL2, SUPavgL2, MEDavgL2, INFavgL1, SUPavgL1, MEDavgL1, IRFdrawsL2, IRFdrawsL1] = VARirbandtrue_endint1shock_historiesEnd2_futureNcons(VAR,OIRF2B_optavg,pick,perc,ndraws,68,typeImpulse,fast,Nconsec);
    [INFavgL2, SUPavgL2, MEDavgL2, INFavgL1, SUPavgL1, MEDavgL1, IRFdrawsL2, IRFdrawsL1] = VARirbandtrue_endint1shock_historiesEnd2_futureNcons(VAR,OIRF2B_optavg,pick,perc,ndraws,68,typeImpulse,1   ,Nconsec);

 
% ---------- keep just the nsteps ahead needed
INFavgL2=INFavgL2(1:nsteps,:,:);
SUPavgL2=SUPavgL2(1:nsteps,:,:);
MEDavgL2=MEDavgL2(1:nsteps,:,:);
INFavgL1=INFavgL1(1:nsteps,:,:);
SUPavgL1=SUPavgL1(1:nsteps,:,:);
MEDavgL1=MEDavgL1(1:nsteps,:,:);
IRFdrawsL2=IRFdrawsL2(1:nsteps,:,:,:);
IRFdrawsL1=IRFdrawsL1(1:nsteps,:,:,:);
% ----------
 pos2plot = find(strcmp(labels,'Industrial production')==1);
%pos2plot = find(strcmp(labels,'SP500 index')==1);
%pos2plot = find(strcmp(labels,'UNC')==1);
MEANavgL2=mean(IRFdrawsL2,4);MEANavgL2=MEANavgL2(1:nsteps,:,:);

FontSize=14;
%% test statistic for the difference between the two state-conditional GIRFs
pctg=68;
select=pick;
%-------------------------------------------------------------------------%
[dIRFinf,dIRFsup,~      ]=VAR_intstat(IRFdrawsL1(1:end-(Nconsec-1),:,:,:),IRFdrawsL2(Nconsec:end,:,:,:),select,pctg);%it computes: IRFdrawsL2 (2 consec shocks)-IRFdrawsL1 (1 shock)

% at the center bands there is the sample difference among responses:
         %LH                  LL (first index refers to state of economy, second index refers to UNCERTANTY)
 dOIRFavg=OIRFavg2B(Nconsec:end,:,pick)-OIRFavg1B(1:end-(Nconsec-1),:,pick);
 
eval(['IRFd_'  num2str(Nconsec) '_low=     dIRFinf(:,pos2plot);']) 
eval(['IRFd_'  num2str(Nconsec) '_hig=     dIRFsup(:,pos2plot);']) 
% Average from bootstrap
 eval(['IRF_'  num2str(Nconsec) '_b  =      MEANavgL2(Nconsec:end,pos2plot,pick);']) 
 eval(['IRFd_' num2str(Nconsec) '_b  =      MEDavgL2( Nconsec:end,pos2plot,pick)  - MEDavgL1( 1:end-(Nconsec-1),pos2plot,pick);'])
 eval(['IRF_'  num2str(Nconsec) '    =     OIRFavg2B( Nconsec:end,pos2plot,pick);']) 
% Difference of point estimates 
 eval(['IRFd_' num2str(Nconsec) '    =     OIRFavg2B( Nconsec:end,pos2plot,pick)  -OIRFavg1B( 1:end-(Nconsec-1),pos2plot,pick)     ;']) 
 if Nconsec == 2
 eval(['IRF_'  num2str(1) '          =                                             OIRFavg1B( 1:end-(Nconsec-1),pos2plot,pick);']) 
% Difference of averages from bootstrap 
                                            MEANavgL1=mean(IRFdrawsL1,4); MEANavgL1=MEANavgL1(1:nsteps,:,:);
 eval(['IRF_'  num2str(1) '_b        =      MEANavgL1(1:end-(Nconsec-1),pos2plot,pick)                                        ;']) 
 end
clearvars -except IRF_* IRFd_* Nconsec

%--------------------%
%StateVar = 'ADS'; 
 StateVar = 'CFNAI'; 
%--------------------% 
 StateType='LevChg';
if     Nconsec == 2
    save(strcat(StateVar,'_',StateType));
elseif Nconsec == 3
    load(strcat(StateVar,'_',StateType));
    save(strcat(StateVar,'_',StateType));
elseif Nconsec == 4
    load(strcat(StateVar,'_',StateType));
    Nconsec=4;
    save(strcat(StateVar,'_',StateType));
end
if Nconsec ==4
%% Figude D.1
FontSize = 46;
COLOR = linspecer(5);
LineWidth =  8;
MarkerSize= 30;

hor = length(IRF_1); zz=zeros(1,hor);
PlotLabel = 'Industrial Production';

figure(111); hold on
%IRF_1 = IRF_1_b;
%IRF_2 = IRF_2_b;
h1=line(1:length(IRF_1), IRF_1, 'color', 'k'       ,'LineStyle','-', 'linewidth',LineWidth);
h2=line(1:length(IRF_2), IRF_2, 'color',COLOR(1,:) ,'LineStyle','--','linewidth',LineWidth,'Marker','s','MarkerSize',MarkerSize);
h3=line(1:length(IRF_3), IRF_3, 'color',COLOR(2,:) ,'LineStyle','--','linewidth',LineWidth,'Marker','o','MarkerSize',MarkerSize);
h4=line(1:length(IRF_4), IRF_4, 'color',COLOR(5,:) ,'LineStyle','--','linewidth',LineWidth,'Marker','d','MarkerSize',MarkerSize);
hold on; plot(1:hor , zz , 'k-' , 'LineWidth' , 2 ); grid; xlim([1 length(IRF_1)]); 
legend([h1 h2 h3 h4],'1 Shock','2 Shocks','3 Shocks','4 Shocks','Location','Best')
ylabel(PlotLabel,'FontSize',FontSize);
xlim([1      hor]);set(gca,'XTick',[6:6:hor],'FontSize',FontSize);
                   set(gca,                  'FontSize',FontSize)
% For paper
  set(gca,'FontSize',52)
  hline = findobj(gcf, 'type', 'line'); set(hline,'LineWidth',6); set(hline(1),'MarkerSize',32); set(hline(3),'MarkerSize',32)
% hAxes = gca;
% reformatTickLabels(hAxes);
  ylim([-0.55 0.15])
%save2pdf('ConsecOutput' ,111,600)


figure(222); hold on
%#CHECK
%plot([IRF_2_b-IRF_1_b,IRF_2-IRF_1])
%plot([IRF_3_b-IRF_2_b,IRF_3-IRF_2])
%plot([IRF_4_b-IRF_3_b,IRF_4-IRF_3])
h2=line(1:length(IRFd_2), IRFd_2, 'color', COLOR(1,:),'LineStyle','--','linewidth',LineWidth,'Marker','s','MarkerSize',MarkerSize);
h3=line(1:length(IRFd_3), IRFd_3, 'color', COLOR(2,:),'LineStyle','--','linewidth',LineWidth,'Marker','o','MarkerSize',MarkerSize);
h4=line(1:length(IRFd_4), IRFd_4, 'color', COLOR(5,:),'LineStyle','--','linewidth',LineWidth,'Marker','d','MarkerSize',MarkerSize);
hold on; plot(1:hor , zz , 'k-' , 'LineWidth' , 2 ); grid; xlim([1 hor]); 
legend([h2 h3 h4],'Incr. due to 2nd shock','Incr. due to 3rd shock','Incr. due to 4th shock','Location','Best')
%ylabel(PlotLabel,'FontSize',FontSize);
xlim([1      hor]);set(gca,'XTick',[6:6:hor],'FontSize',FontSize);
                   set(gca,                  'FontSize',FontSize)
% For paper
  set(gca,'FontSize',52)
  hline = findobj(gcf, 'type', 'line'); set(hline,'LineWidth',6); set(hline(1),'MarkerSize',32); set(hline(2),'MarkerSize',32)
% hAxes = gca;
% reformatTickLabels(hAxes);
 ylim([-.12 0.005]); %set(gca,'YTick',[-0.2:0.1:0 0.05])


% figure(3331); hold on
% hor = length(IRFd_2_hig);
% grpyat=[(1:1:hor)', IRFd_2_low; (hor:-1:1)' IRFd_2_hig(hor:-1:1)]; patch(grpyat(:,1), grpyat(:,2), [0.7 0.7 0.7],'edgecolor', [0.7 0.7 0.7]);
% hold on
% h1= line(1:1:hor, IRFd_2_b, 'color', COLOR(5,:),'LineStyle','--','linewidth',LineWidth);
% hold on; plot(1:hor , zz , 'k-' , 'LineWidth' , 2 ); grid; xlim([1 hor]); 
% xlim([1      hor]);set(gca,'XTick',[6:6:hor],'FontSize',FontSize);
%                    set(gca,                  'FontSize',FontSize)
% % For paper
%   set(gca,'FontSize',52)
%   hline = findobj(gcf, 'type', 'line'); set(hline,'LineWidth',6); set(hline(2),'MarkerSize',32); 
% % hAxes = gca;
% % reformatTickLabels(hAxes);
% ylim([-.2 0.01]); set(gca,'YTick',[-0.2:0.1:0 0.05]) 
% 
% figure(3332); hold on
% hor = length(IRFd_3_hig);
% grpyat=[(1:1:hor)', IRFd_3_low; (hor:-1:1)' IRFd_3_hig(hor:-1:1)]; patch(grpyat(:,1), grpyat(:,2), [0.7 0.7 0.7],'edgecolor', [0.7 0.7 0.7]);
% hold on
% h1= line(1:1:hor, IRFd_3_b, 'color', COLOR(5,:),'LineStyle','--','linewidth',LineWidth);
% hold on; plot(1:hor , zz(1:hor) , 'k-' , 'LineWidth' , 2 ); grid; xlim([1 hor]); 
% xlim([1      hor]);set(gca,'XTick',[6:6:hor],'FontSize',FontSize);
%                    set(gca,                  'FontSize',FontSize)
% % For paper
%   set(gca,'FontSize',52)
%   hline = findobj(gcf, 'type', 'line'); set(hline,'LineWidth',6); set(hline(2),'MarkerSize',32); 
% ylim([-.2 0.01]); set(gca,'YTick',[-0.2:0.1:0 0.05]) 
%  
% figure(3333); hold on
% hor = length(IRFd_4_hig);
% grpyat=[(1:1:hor)', IRFd_4_low; (hor:-1:1)' IRFd_4_hig(hor:-1:1)]; patch(grpyat(:,1), grpyat(:,2), [0.7 0.7 0.7],'edgecolor', [0.7 0.7 0.7]);
% hold on
% h1= line(1:1:hor, IRFd_4_b, 'color', COLOR(5,:),'LineStyle','--','linewidth',LineWidth);
% hold on; plot(1:hor , zz(1:hor) , 'k-' , 'LineWidth' , 2 ); grid; xlim([1 hor]); 
% xlim([1      hor]);set(gca,'XTick',[6:6:hor],'FontSize',FontSize);
%                    set(gca,                  'FontSize',FontSize)
% % For paper
%   set(gca,'FontSize',52)
%   hline = findobj(gcf, 'type', 'line'); set(hline,'LineWidth',6); set(hline(2),'MarkerSize',32); 
% ylim([-.2 0.01]); set(gca,'YTick',[-0.2:0.1:0 0.05])   
end